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Bistability induced by nonlinear Kerr effect in arrays of coupled waveguides is studied and shown 
to be a means to conceive light detectors that switch under excitation by a weak signal. The detector 
is obtained by coupling two single ID waveguide to an array of coupled waveguides with adjusted 
indices and coupling. The process is understood by analytical description in the conservative and 
continuous case and illustrated by numerical simulations of the model with attenuation. 
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Introduction. Arrays of adjacent optical dielectric 
waveguides coupled by power exchange between guided 
modes p] Q , have allowed to conceive devices possessing 
extremely rich optical properties, already at the linear 
level, such as anomalous refraction and diffraction 

At nonlinear level, for intensity-dependent refractive 
index (optical Kerr effect), these waveguide arrays be- 
come soliton generators EL as experimentally demon- 
strated in |ol la. 171 Isl M. llflTll| . The model is a discrete 
nonlinear Schrodinger equation (NLS) and nonlinearity 
then manifests by self-modulation of an input signal (in- 
jected radiation) that propagates as a NLS discrete soli- 
ton |j| . These systems possess also intrinsically dis- 
crete properties |l3| a nd the geometry can be varied to 
manage dispersion 14| . 

A fundamental property of nonlinear systems that has 
not been considered in waveguides arrays is the bista- 
bility induced by nonlinearity. It is the purpose of this 
work to propose and study a device where bistability, and 
consequent switching properties, could be observed and 
used to conceive for instance a detector of light sensitive 
to extremely weak signal. 

We shall make use of the possibility to drive a waveg- 
uide array, in the forbidden band gap, through directional 
coupling by boundary waveguides, which results in the 
generation of (discrete) gap solitons [lfij . produced by 
nonlinear supratransmission |17|. and used to conceive 
resonators with nonlinear eigenstates • 

The waveguide array of figure ^ m nonlinear Kerr 
regime, driven by two single waveguides of index no > n\, 




is operated in such a way that the guided modes in the 
two lateral waveguides would not linearly propagate in the 
array due to index difference. Nonlinear induced bista- 
bility allows then to adjust the input radiation such as 
to work at the bifurcation point and switch from a state 
of vanishing intensity output (from the central waveg- 
uides) to a state of strong intensity output, as displayed 
on fig|2]that represents the intensity of the flux in the ar- 
ray when operated close to the bifurcation point. As the 
switch can be operated by a very weak signal (0.03% of 
the input amplitude for fig|2Jl, the device is a candidate 
for an ultrasensitive bistable detector. 
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FIG. 1: Waveguide array (index profile) driven by 2 lateral 
waveguides of index no > ni . The arrows symbolize the in- 
jected radiation. 
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FIG. 2: (Color inline) Flux intensity along the waveguides 
driven at the threshold of bistability. This is the result of 
numerical simulations of the model JSJ with N — 4, Q = 
Qo = 4, 7 = 70 = 0.001, V = 0.5, and input amplitudes 
b — 3.0946 (large arrows). In the lower plot the small arrow 
represents a signal of amplitude 0.001 enough to trigger the 
switch. 

After having briefly recalled the derivation of the dis- 
crete NLS model for the envelope of the guided modes 
and the proper boundary conditions that results from the 
physical context of figE we describe theoretically and 
numerically the bistability process. The theory calls to 
the continuous limit which gives an accurate description 
of the discrete case with enough waveguides and strong 
coupling. 

The model. Assuming instantaneous response of the 
Kerr medium (of nonlinear susceptibility x) an( i pl &_ 
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nar wave guides along the vertical direction x, the elec- 
tric field can be sought as E{r,t) — u x (E(y,z) e~ tuJt + 
c.c.)/y/3x and Maxwell's equations reduce to 



d 2 



d 2 , „ lo 2 



dy 2 dz' 



) E+^(n 2 +\E\ 2 )E = 0, 



(1) 



where the linear index n — n(r) may vary in space. Inter- 
ested in the stationary regime, the fast oscillating nonlin- 
ear term has been discarded (assuming no phase match- 
ing between third harmonics). 

We consider the index profile of fig^Jand assume a field 
component E(y, z) in the form 



E{y,z) = e^2,£ j {ey)'4> j {e 2 z)e 



i/3z 



(2) 



where £j is a guided mode, ipj its amplitude and e a small 
parameter to be defined later. The slow variation along z 
results from the coupling between adjacent waveguides. 
Guided modes in both lateral waveguides, that do not 
propagate transversaly in the array, result from the nec- 
essary condition 



f] 2 - -^n 2 < 0, /3 2 - -^nj > 0. 



(3) 



Then, following > we insert J2J) in , integrate over y 
and obtain the coefficients that couple the driving lat- 
eral wave guides to the array (Co), the array waveg- 
uides together (Ci), and that result from the overlap- 
ping of modes £j across the separation between waveg- 
uides. Upon defining the dimensionless variable z' — 
zuj 2 I {2[3c 2 ) the coupling coefficients are scaled by 



7C1, Qo 



c 2 Cl 



and we define the essential constant 



V = n 
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2Q, V>0, V ~ e 2 



(4) 



(5) 



This defines the small parameter e used in expansion J5| , 
and provides the necessary and sufficient condition (V > 
0) for evanaescent waves along transverse direction in the 
array. 

We obtain eventually the model (forgetting the prime 
on z and setting e = 1 to go back to physical units) 



i(d z + 70)^0 = Qoipi, i(d z + jo)ipN = Qo4>n-i 
i{d z + 7)V> 3 + Q (-0J-1 + i>j+i ~ 2-0j) + l^'l 2 ^ = 



(6) 



where j runs from j = 1 to j = N — 1. Attenuation has 
been included (imaginary parts 7 and 70 of the dielectric 
constants) to account for a realistic physical situation, 
and the actual amplitudes in the lateral driving waveg- 
uide have been rescaled by replacing Co^o.jv by Citpo t N- 
Prescription of the injected energy flux in the lateral 
waveguides results in defining the amplitudes ipo and tpN 
in z = (say b) , which are associated to a set of vanishing 



ipj (0) for the rest of the array. The initial condition is 
thus 

V>o(0) =iM0) = b, ^(0) = 0, j = l---N-l (7) 

which complete equations JJjJ. The point is that non- 
linearity causes the existence of a threshold b s allowing 
to take advantage of bistabilility of stationary states, as 
described now analytically and numerically. 

Analytic description. The continuous limit of © is 
obtained for 



y 



(8) 



(note that this dimension y is not the physical dimension 
y of expression |2"jl) and reads (without attenuation) 



ilp z +1pyy + \lp\ 2 1p = Vlp. 



(9) 



The boundary data that represent the injected flux are 
taken as constants, namely 



^(0,*)=V(£,«) = &, 



(10) 



As far as the light tunnels from the driving waveguides 
to the array, the amplitude of light in the driving waveg- 
uides decays along z. However, if the coupling Qo and 
the attenuation 70 are small enough, the intensities in the 
driving waveguides decay slowly along z and the above 
boundary-value problem quite correctly matches numer- 
ical simulation of the discrete model iJBJ . 

The stationary solutions of the above problem are ob- 
tained by assuming a real valued solution depending on 
y only, namely ip(y, z) — ip(y) solution of 



(11) 



From the symmetry of the above bounds, a uniformly 
bounded function requires that there exists yo such that 



^3/(2/0) = 0, ijj{y a )=B. 



(12) 



The values b and B are referred respectively to as input 
and output amplitudes of the field. By integrating l|ll|) 
with l|12l) . one gets 



Wv) 2 



-(B 2 
2 V 



^ 2 ){^ 2 + B 2 -2V). 



(13) 



This equation has solutions of different types |20| de- 
pending on the relative values of B and V . In particu- 
lar, an input value 6 of the boundary driving can pro- 
duce different output amplitudes B. This property holds 
when the input amplitude b is less than a value b s , called 
the supratransmission threshold, above which the bound- 
ary driving induces an instability which generates soliton 
emission in the array |2l| . 

Case I: B 2 > 2V. The integration of lO on [0, L] 
produces the relation 



6 = B-cn[/3L/2,k], 



(14) 
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where cn[-,k] is the cosine-amplitude Jacobi elliptic func- 
tion of modulus k. The parameters j3 and k obey the 
relation 



(3 2 = B 2 - V, k 2 



D 2 



(15) 



2{B 2 -Vy 

while the whole solution 1)11(1 has the following form 

^(2/) = B-cn[/3(y-L/2),k]. (16) 

Case II: V < B 2 < 2V. The integration of JUJ gives 
here 



V) 



6 = S-dn[/3L/2,k], 
^ // (y) = S-dn[/3(y-L/2),k] 



Case III: B 2 < V . In this case we obtain 



b = 



B 



dn[/3L/2,k]' 
IV - B 2 



k 2 = 
B 



2{V-B 2 ) 
2V- B 2 



dn[/3(y-L/2),k]- 



(17) 

(18) 
(19) 

(20) 

(21) 
(22) 



This last solution has a threshold amplitude b = b s 
reached when the function b(B) in (|20|l has a vanishing 
derivative. The related expression of b s is not explicit al- 
though easily evaluated numerically. Note that in general 
one has a finite set of thresholds (their number depend 
on the lenght L) but only the first one is of interest for 
our purpose. Note also that the above three solutions 
reach their maximum amplitude B in the array center 

y = L/2. 

Bistable behavior. The numerical simulations of sys- 
tem ©, now compared to the above analytical solu- 
tions, are performed with N = 18, Q = 20, Qq = 0.05 
and V = 0.5. The system locks to a stationary so- 
lution thanks to the attenuation factors 7 = 0.02 and 
70 = 0.002. In order to avoid an initial shock, the cou- 
pling Q is actually smoothly set by using Qo(l — e _z / 5 ). 
The input-output dependence resulting from of (|14fl , i|17|) 
and H20fl is plotted on fig[3]and the profiles of the differ- 
ent exact solutions ifTCI) . (|T§|) and lt22*|) corresponding to 
the different output amplitudes are presented in fig0] 

The regimes of stationary solutions are reached by nu- 
merical simulations of JfJJl as follows. First we inject into 
the driving waveguides the beams with the initial ampli- 
tude for b < b s . This produces the type III solution, 
a nonlinear analogue to the evanescent wave profile. As 
the amplitudes l^o^)! and |^>Ar(z)| in the driving waveg- 
uides decay along z we monitor the amplitude \ipN/2{z)\ 
at the middle of the array as a function of the driv- 
ing amplitudes. At some distance zq from the origin, 
the system stabilizes to a stationary profile (as shown 
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FIG. 3: Hysteresis loop for the input-output amplitude values 
obtained from the analytical solutions of type I Ijl4[l . 77 l|17[l 
and III (I20t respectively. Crosses are the results of numerical 
simulations. The dashed vertical line represents the value of 
the supratransmission threshold b s , and ther two dashed lines 
set upper and lower limits for the appearance of type II solu- 
tion. Crossing of the full vertical line with the hysteretic curve 
defines input-output values 1, 2 and 3 chosen for /tgO[ 
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FIG. 4: Analytic solutions (dashed lines) corresponding to 
the single input amplitude b = 0.18 compared to numerical 
experiments (solid lines). 



by figEJ) and we measure there the corresponding ampli- 
tudes B = \"4>n/2{ z o)\ and b = |^o(zo)|- plotted then as 
the crosses along the lower branch (III) in figOJ (Note 
that only the solutions of type I and III can be reached 
by numerical simulations, while the type II is unstable 
and eventually decays to the stable solution III. Note 
also that longer lengths would allow for multiple output 
amplitudes, which we do not consider here.) 

The same procedure is used when injecting the beams 
with b > b s . After a regime of gap solitons emission, the 
amplitudes along the lateral waveguides decrease below 
the threshold b s , allowing the system to lock to the solu- 
tion of type I. The obtained numerical values are plotted 
(crosses) along the upper branch (I) of the hysteresis loop 
in figGl 

The analytical curves slightly overestimate the numer- 
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FIG. 5: Evolution of input (inset) and output (main plot) am- 
plitudes along the propagation direction z. Note the difference 
between output amplitudes \ip(L/2, z)\ in two regimes although 
the initial values of the amplitudes in the boundary waveguide 
b slightly differ from each other. 

ical results because of the attenuation (optical losses in 
the laboratory experiments) included in our numerical 
simulations. Another consequence of the attenuation 
is the existence of a second (nonzero) threshold input 
amplitude where the solution bifurcates back from the 
regime I to the regime III. For a vanishing attenuation, 
this second threshold would be exactly zero according to 
the analytical solution (|l(j[l . The existence of this lower 
end threshold is clearly of importance in view of experi- 
mental realization but, as being directly related to atten- 
uation, it depends on the precise physical context. Let us 
also remark that the output amplitude is sensitive to the 
input amplitude values (in both regimes). In particular 
the curves (and dots) of figEJshow that this dependence 



is roughly linear as the slopes of the functions B(b) in 
regimes I and III have approximately the same slopes. 

As a matter of fact we display on the figElthe evolution 
of the amplitudes l^oO^OI = \iPn(z)\ along the driving 
waveguides and the amplitude in the middle of the ar- 
ray \'4 > n/2(z) \ obtained from a numerical simulation. The 
difference between the two output signals at y = L/2 in 
these two regimes differs substantially (approximately by 
two orders of magnitude). 

The existence of bistable regimes that switch from a 
low to a high output is a candidate for an ultra sensi- 
tive bistable light detector. Moreover such a property 
holds also in the fully discrete case as displayed in fig|3 
(though we do not have in that case an analytical descrip- 
tion). Particularly, if one injects in the lateral waveguides 
beams of amplitude slightly below the supratransmission 
threshold b s , a regime of type III establishes (with small 
output amplitude). Then any perturbation of the array 
by means of a weak input signal destabilizes the system 
which, after transient emission of gap solitons, reaches a 
type I regime with large output amplitude. This is pre- 
cisely the process described in the intensity plot of fig[2l 
where the signal producing the switch carries 0.03% of 
the driving amplitude. 

Note that in such a fully discrete case, the device works 
as a digital detector, triggered by any signal received by 
the central waveguides. However in the continuous case, 
the output signal amplitude can be (approximately) lin- 
early related to the input amplitude and, if the variable 
parameter is the input amplitudes in the lateral waveg- 
uides, then the device can be operated as an amplificator. 
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